# Example input file to test the neo-hookean material
# (initially based on the tensor_mechanics tutorial, Step 1,
# part 2: 2D simulation of uniaxial tension with linear
# elasticity; Step 2, part2, and
# moose/modules/tensor_mechanics/test/tests/
# finite_strain_elastic/finite_strain_elastic_new_test.i)

# Status: same values as FEBio

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  order = first
[]

[Mesh]
  file = cube.msh
  dim = 3
[]

[Modules/TensorMechanics/Master]
  [./all]
    # Finite strain is useful for hyperelasticity
    strain = FINITE
    # Detects the change of the mesh to second order and
    # automatically sets the variables
    add_variables = true
    # Automatically creates the auxvariables and auxkernels
    # needed to output these stress quanities
    generate_output = 'strain_xx strain_xy strain_xz strain_yx strain_yy strain_yz strain_zx strain_zy strain_zz stress_xx stress_xy stress_xz stress_yx stress_yy stress_yz stress_zx stress_zy stress_zz'
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeElasticityTensorNeohookean
    # # E = 30e6, nu = 0.499
    # lambda = 4.993328885923944E+9
    # mu = 1.000667111407605E+7
    E = 30e6
    nu = 0.3
    # # To check the Voigt notation
    # lambda = 2
    # mu = 5
  [../]

  [./stress]
    type = ComputeStressNeohookean
    # # E = 30e6, nu = 0.499
    # lambda = 4.993328885923944E+9
    # mu = 1.000667111407605E+7
    E = 30e6
    nu = 0.3
    # # To check the Voigt notation
    # lambda = 2
    # mu = 5
  [../]
[]

[BCs]
  [./fix_x]
    type = PresetBC
    variable = disp_x
    boundary = left
    value = 0
  [../]
  [./fix_y]
    type = PresetBC
    variable = disp_y
    boundary = back
    value = 0
  [../]
  [./fix_z]
    type = PresetBC
    variable = disp_z
    boundary = bottom
    value = 0
  [../]
  [./pull]
    # Make displacement equal to time
    function = 't'
    type = FunctionPresetBC
    variable = disp_z
    boundary = top
  [../]
  # [./pull]
  #   type = PresetBC
  #   variable = disp_z
  #   boundary = front
  #   value = 10e-6
  # [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

# [Executioner]
#   type = Transient
#   end_time = 500e-6             # = \Delta{l}, stretch = 1.5
#   num_steps = 10
#
#   solve_type = 'NEWTON'
#
#   petsc_options = '-snes_ksp_ew -ksp_view_mat'
#   petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap -ksp_gmres_restart'
#   petsc_options_value = 'asm lu 1 101'
# []
[Executioner]
  type = Steady

  solve_type = 'NEWTON'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap -ksp_gmres_restart'
  petsc_options_value = 'asm lu 1 101'
[]

[AuxKernels]
  [./C1111]
    type = RankFourAux
    rank_four_tensor = elasticity_tensor
    index_i = 0
    index_j = 0
    index_k = 0
    index_l = 0
    variable = C1111_aux
  [../]
  [./C1122]
    type = RankFourAux
    rank_four_tensor = elasticity_tensor
    index_i = 0
    index_j = 0
    index_k = 1
    index_l = 1
    variable = C1122_aux
  [../]
  [./C2323]
    type = RankFourAux
    rank_four_tensor = elasticity_tensor
    index_i = 1
    index_j = 2
    index_k = 1
    index_l = 2
    variable = C2323_aux
  [../]
  [./J]
    # Loosely based on https://groups.google.com/d/msg/moose-users/CLGXpbQu7Ko/W0HwiWoCDAAJ
    type = RankTwoScalarAux
    execute_on = timestep_end
    # Take a tensor which is defined in your kernels (source code)
    rank_two_tensor = deformation_gradient
    # Give it a name (that you define in the AuxVariables)
    variable = detF
    # Get the determinant (third invariant) (For more, look
    # at [moose/modules/tensor_mechanics/src/auxkernels/
    # RankTwoScalarAux.C] and [moose/modules/tensor_mechanics/
    # src/utils/RankTwoScalarTools.C])
    scalar_type = ThirdInvariant
  [../]
[]

[AuxVariables]
  # Main components of the elasticity tensor (all other
  # values are zero)
  [./C1111_aux]
    # Value of the diagonal elements of the upper-left 3x3
    # sub-matrix
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./C1122_aux]
    # Value of the off-diagonal elements of the upper-left
    # 3x3 sub-matrix
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./C2323_aux]
    # Value of the diagonal elements of the lower-right
    # 3x3 sub-matrix
    order = CONSTANT
    family = MONOMIAL
  [../]

  # Determinant of the deformation gradient tensor
  [./detF]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Postprocessors]
  # See other postprocessors here:
  # [moose/framework/src/postprocessors/]
  [./ave_stress_bottom_z]
    type = SideAverageValue
    variable = stress_zz
    boundary = bottom
  [../]
  [./ave_strain_top_z]
    type = SideAverageValue
    variable = strain_zz
    boundary = top
  [../]
  [./el_stress_zz]   # This could also be stress_zz
    # This post-processor does sum over element and divides
    # by volume. When using this post-processor, it is
    # better to set allow_renumbering = false in the [Mesh]
    # block
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  [../]
  [./el_stress_yy]
    # This post-processor does integration over the element,
    # and normalises by volume (and does not need elementid)
    type = ElementAverageValue
    variable = stress_yy
  [../]
  [./ave_C1111]
    type = ElementalVariableValue
    variable = C1111_aux
    elementid = 0
  [../]
  [./ave_C1122]
    type = ElementAverageValue
    variable = C1122_aux
  [../]
  [./ave_C2323]
    type = ElementalVariableValue
    variable = C2323_aux
    elementid = 0
  [../]
  [./det_F]
    type = ElementAverageValue
    variable = detF
  [../]
  [./disp_top_z]
    type = SideAverageValue
    variable = disp_z
    boundary = top
    execute_on = 'initial timestep_end'
  [../]
  [./cube_stretch]
    # Use a MOOSE function to calculate the coupon of the
    # stretch
    type = FunctionValuePostprocessor
    # This function must exist in the [Functions] block
    function = overal_stretch_z
    execute_on = 'initial timestep_end'
  [../]
[]

[Functions]
  # Create a function to calculate the instantaneous stretch
  # as measured from the top to the bottom of the coupon
  # (remember that this is 1/8 of a coupon)
  [./overal_stretch_z]
    type = ParsedFunction
    # The variables and its values are given in something
    # similar to a table (we take disp_top_z from the
    # [Postprocessors] block; the extra space is for
    # readability):
    vars = "l0    displ"
    vals = "1e-3  disp_top_z"
    # This is the value that the ParsedFunction gets (both
    # variables are defined within this block)
    # \(Stretch = \Delta{l} / l_{0}\)
    value = "(displ) / l0 + 1"
  [../]
[]

[Outputs]
  [./o]
    type = Exodus
    # May be I can get values at the nodes with ParaView in
    # this way?
    elemental_as_nodal = true
  [../]
  # exodus = true
  print_perf_log = true
  # Do not print data on screen
  console = false
  # Export CSV (comma-separated data as text file)
  csv = true
[]
